/***************************************************************************
 *
 *   Copyright (C) 2003 by Haydon Knight
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/

#include "environ.h"

#include "dsp/Transformation.h"
#include "dsp/TimeSeries.h"
#include "dsp/BitSeries.h"

#include "dsp/ChannelOrder.h"

using namespace std;

dsp::ChannelOrder::ChannelOrder()
  : Transformation<BitSeries,TimeSeries>("ChannelOrder",outofplace)
{
  rapid = Polarisation;
}

dsp::ChannelOrder::~ChannelOrder(){ }

void dsp::ChannelOrder::transformation()
{
  if( verbose )
    cerr << "dsp::ChannelOrder::transformation with input->ndat=" 
         << input->get_ndat() << endl;

  if(input->get_nbit()!=32)
    throw Error(InvalidState,"dsp::ChannelOrder::transformation()",
		"nbit!=32- this routine is only set up for floats generated by dsp::TimeOrder\n");

  output->Observation::operator=( *input );
  output->set_swap( false );  
  output->resize( input->get_ndat() );

  if( rapid==Channel ){
    // number of floats between (t0,f0) and (t1,f0) of a BitSeries
    const unsigned input_stride = input->get_nchan()*input->get_npol()*input->get_ndim();
    const unsigned output_stride = output->get_ndim();
    
    for( unsigned idim=0; idim<output->get_ndim(); idim++){
      for( unsigned ichan=0; ichan<output->get_nchan(); ichan++){
	for( unsigned ipol=0; ipol<output->get_npol(); ipol++){
	  float* in = (float*)input->get_rawptr() + idim + input->get_ndim()*         (ichan + input->get_nchan()*ipol); 
	  float* out = (float*)output->get_datptr(ichan,ipol)+idim;
	  
	  const unsigned isamp_end = output->get_ndat()*output->get_ndim(); 
	  unsigned input_samp=0; 
	  
	  for( unsigned isamp=0; isamp<isamp_end;
	       isamp+=output_stride, input_samp+=input_stride)
	    out[isamp] = in[input_samp];
	}
      }
    }    

  }
  else if( rapid==Polarisation ){
    // number of floats between (t0,f0) and (t1,f0) of a BitSeries
    const unsigned input_stride = input->get_nchan()*input->get_npol()*input->get_ndim();
    const unsigned output_stride = output->get_ndim();
    
    for( unsigned idim=0; idim<output->get_ndim(); idim++){
      for( unsigned ipol=0; ipol<output->get_npol(); ipol++){
	for( unsigned ichan=0; ichan<output->get_nchan(); ichan++){
	  float* in = (float*)input->get_rawptr() + idim + input->get_ndim() * (ipol + input->get_npol()*ichan);
	  float* out = (float*)output->get_datptr(ichan,ipol)+idim;
	  
	  const unsigned isamp_end = output->get_ndat()*output->get_ndim(); 
	  unsigned input_samp=0;

	  for( unsigned isamp=0; isamp<isamp_end;
	       isamp+=output_stride, input_samp+=input_stride)
	    out[isamp] = in[input_samp];
	}
      }
    }

    if( verbose )
      fprintf(stderr,"In ChannelOrder::transformation() rapid=pol\n");

  }

  if( verbose )
    cerr << "dsp::ChannelOrder::transformation exit with output->ndat="
         << output->get_ndat() << endl;
}

